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Abstract 


This  report  describes  the  background,  theory  and  simulation  of  the  internal  wave  wake 
from  ships.  Realistic  models  of  both  internal  layers  and  hulls  are  employed.  The  internal 
layer  models  are  based  on  published  vertical  profiles  of  the  Brunt-Vaisala  frequency. 
Two  basic  profiles  are  employed.  The  first  is  a  profile  appropriate  to  the  Strait  of 
Gibraltar  and  is  modeled  using  a  Lorentzian  distribution;  the  parameters  are  the  peak 
Brunt-Vaisala  frequency,  its  depth  and  the  width  of  the  distribution.  If  required,  these  can 
be  altered  to  apply  to  other  situations.  The  second  profile  is  appropriate  to  that  existing 
during  the  Loch  Linnhe  trials  in  Scotland. 

The  hull  models  represent  warship  and  merchant  ship  hulls.  The  warship  hull  model  is  a 
modified  Taylor  model  and  the  block  coefficient  can  be  adjusted  continuously  within 
suitable  limits;  the  parent  model  was  based  on  a  twin  screwed  cruiser.  The  merchant  ship 
model  was  developed  by  the  David  Taylor  Model  Basin  and  applies  to  single  screwed 
ships.  There  are  5  separate  models  covering  a  range  of  block  coefficients.  The  Wigley 
model,  which  resembles  a  canoe,  a  single  source  model  and  a  simple  source  and  sink 
model  (doublet)  are  available  in  part  for  diagnostics. 

The  simulations  successfully  reproduce  some  observations  of  wakes  during  the  Loch 
Linnhe  trials,  in  which  wakes  with  significant  structure  were  generated.  Wakes  generated 
on  the  Gibraltar  layer  appear  to  be  on  the  edge  of  detectability  by  airborne  or  space-bome 
high  resolution  radar  and  consist  mainly  of  a  V-shape  with  little  structure  in  between  the 
arms. 


Executive  Summary 


A  full  understanding  of  internal  wave  wakes  is  required  if  these  wakes  are  to  be  useful  for 
maritime  surveillance.  The  objective  of  this  work  is  to  gain  this  understanding  by 
numerically  simulating  wakes  for  realistic  internal  layers  using  realistic  hull  models 
moving  at  operational  speeds.  Simplified  layer  models  can  be  useful  and  previous  work 
has  shown  that  rough  agreement  with  observations  can  be  achieved  under  some 
circumstances.  Realistic  models  offer  the  opportunity  for  more  reliable  conclusions  about 
wake  visibility,  for  example  in  radar  imagery.  However,  the  present  study  is  limited  to  the 
super-critical  case  when  the  ship  speed  is  greater  than  the  maximum  wave  velocity. 

Previous  layer  models  have  included  the  discrete  type  where  a  fluid  of  low  density,  such 
as  fresh  water,  lies  on  a  fluid  of  high  density,  such  as  salt  water,  forming  a  discrete 
horizontal  interface  on  which  waves  can  propagate.  At  the  next  level,  the  layer  can 
comprise  a  thick  layer  of  constant  Brunt-Vaisala  (B-V)  frequency  sandwiched  between 
fluids  of  constant  but  differing  densities.  In  practice  the  B-V  frequency  changes 
continuously  and  may  exhibit  several  peaks.  The  presence  of  irregularities  in  the  B-V 
frequency  and  its  slow  fall-off  at  depth  introduce  some  technical  problems  in  the 
numerical  calculations.  Solutions  to  these  problems  are  addressed. 

The  simulations  have  been  validated  by  comparing  the  results  with  observations  from  the 
Loch  Linnhe  trials.  Agreement  with  the  simulated  and  observed  crest  patterns  is  excellent 
when  the  published  Loch  Linnhe  vertical  profile  is  used.  For  the  ships  in  the  trials,  the 
predicted  amplitudes  are  also  of  the  correct  magnitude;  unlike  for  the  discrete  layer 
model,  no  adjustments  are  needed. 

Naturally  occurring  internal  waves  in  the  ocean  are  very  common  in  littoral  areas  and  can 
be  observed  using  satellite  borne  optical  sensors  and  high  resolution  radar.  Their 
existence  implies  the  presence  of  significant  internal  layers.  The  layer  characteristics  are 
described  in  the  published  literature.  In  the  open  ocean,  natural  internal  waves  are  less 
common  and  there  is  a  dearth  of  information.  However,  layers  should  exist,  especially 
during  summer  in  the  tropics. 

The  wake  and  its  detectability  depend  on  a  number  of  factors.  The  principal  factors 
associated  with  the  ship  appear  to  be  its  principal  dimensions  (length,  beam  and  draft),  its 
block  coefficient  and  its  speed.  The  principal  factors  associated  with  the  layer  are  its 
maximum  B-V  frequency,  the  depth  and  width  of  the  layer  and  the  presence  and  width  of 
a  mixed  layer  at  the  surface.  The  latter  tends  to  suppress  a  wake  and  the  presence  of 
mixed  layers  probably  explains  why  internal  wave  wakes  are  not  frequently  observed  in 
radar  imagery. 

For  layers  appropriate  to  littoral  ocean  conditions,  the  results  suggest  that  internal  wave 
wakes  generated  by  shipping  should  commonly  occur  but  are  usually  on  or  below  the 
threshold  of  detectability  using  present  radar  technology.  There  is  some  expectation  that 


conditions  can  be  identified  when  wakes  can  be  detected  reliably  and  that  improvements 
to  radar  and  other  systems  might  be  made  to  achieve  a  useful  probability  of  detection. 
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Introduction 

As  noted  in  previous  reports  [1]  [2],  ship  wakes  and  particularly  those  in  space-borne 
high  resolution  radar  images  can  be  important  for  maritime  surveillance.  They  can  yield 
information  about  the  ship  course  and  speed  and  help  to  verify  the  data  from  other 
sensors,  such  as  Automatic  Identification  Systems  (AIS),  which  must  be  carried  by  large 
ships  according  to  international  maritime  law.  The  purpose  of  this  report  is  to  examine 
the  surface  manifestation  of  internal  wave  wakes  generated  by  ships  using  realistic 
models  of  both  internal  layers  and  ship  hulls  through  numerical  simulation.  This  is  in 
contrast  to  previous  reports  and  papers  for  which  the  model  of  the  undisturbed  layer  was 
either  a  discrete  horizontal  planar  interface  bounding  fluids  of  different  density  [3]  or 
involved  a  few  sections  in  which  the  density  varied  linearly  with  depth  (e.g.  [4]).  In  some 
cases  the  focus  has  been  on  the  crest  pattern  (e.g.  [6]). 

This  study  is  limited  to  the  super-critical  case  when  the  ship  speed  is  greater  than  the 
maximum  speed  of  internal  waves,  which  is  typically  less  than  1  m/s  and  often  much  less. 
In  the  super-critical  regime,  which  corresponds  to  operational  speeds,  the  wake  comprises 
only  divergent  waves  with  no  transverse  wave  structure.  As  the  ship  speed  increases,  the 
wake  becomes  narrower.  At  speeds  lower  than  the  maximum  wave  speed,  the  wake 
comprises  both  divergent  and  transverse  waves  and  its  crest  pattern  resembles  that  of  the 
Kelvin  wake  but  the  wake  angle  is  greater  than  19.5°.  The  presence  of  transverse  waves 
can  create  a  large  drag  on  the  ship  and  this  is  known  as  the  “dead  water”  effect.  Thus  it  is 
possible  that  the  low  speed  case  could  also  be  important  for  wake  detection. 

The  Loch  Linnhe  trials  were  conducted  between  1989  and  1994  and  provided 
documented  internal  wave  wakes  including  radar  images.  Fortunately  comprehensive 
documentation  for  one  or  two  wakes  is  available  in  the  open  literature  and,  as  in  [2],  this 
forms  the  observational  data  for  comparison  with  simulations.  Simulation  of  the  high 
resolution  radar  image  of  a  ship  wake  located  on  the  surface  of  the  sea  requires  an 
understanding  of  hydrodynamics,  radar  scattering  and,  for  space-bome  applications, 
Synthetic  Aperture  Radar  (SAR)  processing.  To  simplify  the  process,  only  the 
hydrodynamics  is  considered  here.  Radar  imagery  of  gravity  wave  wakes  (this  class 
includes  internal  wave  wakes)  is  affected  by  horizontal  surface  fluid  flows  associated 
with  wakes  and  these  are  compared  with  those  derived  from  measurements  taken  during 
the  trials. 

During  the  trials  the  sub-surface  fluid  flows  were  monitored  as  a  function  of  time  at 
various  depths  on  a  vertical  line  at  a  location  near  the  ship’s  track  [5].  In  [2],  using  the 
discrete  layer  model  it  was  shown  that  the  crest  pattern  roughly  matched  the  observed 
pattern  and  that  surface  flows  of  approximately  the  correct  magnitude  were  predicted. 
However,  it  was  not  possible  to  match  both  the  flow  velocities  and  the  number  of 
oscillations  in  the  time  record  simultaneously.  The  number  of  oscillations  could  be 
matched  by  increasing  the  assumed  depth  of  the  discrete  interface  but  then  the  amplitudes 
were  significantly  reduced.  As  will  be  seen,  this  problem  is  overcome  when  the 
appropriate  diffuse  layer  is  modeled. 
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The  theory  and  implementation  of  the  simulation  of  internal  wave  ship  wakes  has  been 
sketched  in  the  literature  [4]  and  indeed  a  correspondence  between  simulations  and 
internal  wave  wakes  has  been  claimed,  though  a  simplified  layer  model  was  employed. 
Realistic  internal  layer  models  are  based  on  the  Brunt- Vaisala  (B-V)  frequency  vertical 
profile.  The  B-V  frequency  is  the  natural  frequency  of  oscillation  of  a  parcel  of  fluid  that 
is  displaced  by  a  small  vertical  distance.  In  [4]  the  B-V  frequency  was  constant  over  the 
layer  and  zero  outside  it  so  that  a  rectangular  vertical  profile  was  used.  One  purpose  of 
this  report  is  to  describe  the  simulation  process  in  detail  (using  a  realistic  profile)  and  to 
examine  the  effects  of  varying  the  parameters  of  the  layer,  the  ship  hull  shape  and  the 
ship  speed  on  the  appearance  of  the  wake  and  on  the  magnitude  of  the  surface  flows. 

Large  amplitude  internal  waves  can  develop  into  internal  shocks  and  solitons;  the  waves 
can  propagate  over  long  distances  while  their  shape  hardly  changes.  Larger  waves  may 
break.  These  are  non-linear  effects.  For  small  amplitude  waves  a  linear  theory  may 
suffice  and  this  should  usually  be  appropriate  for  ship  generated  waves  especially  in  the 
far  wake.  The  correspondence  between  observed  internal  wake  crest  patterns  and  the 
theoretical  predictions  suggests  that  a  linear  theory  is  adequate.  If  ships  were  to  generate 
soliton  wakes,  their  wakes  would  be  primarily  in  the  fonn  of  a  ‘V’  with  little  internal 
structure. 

Wake  formation  can  be  understood  in  terms  of  gravity  waves  emanating  from 
disturbances  created  by  the  movement  of  the  hull  through  the  water  as  it  first  pushes  it 
aside  and  as  it  subsequently  flows  back  into  its  previous  position.  If  turbulence  is 
neglected,  the  disturbances  are  stationary  in  the  reference  frame  of  the  ship  and  this 
considerably  simplifies  the  theory.  Another  simplification  is  possible  if  we  focus  on  the 
far  wake  because  the  waves  are  dominated  by  dispersive  effects  and  a  wave  packet  that 
travels  a  large  distance  becomes  closely  sinusoidal,  regardless  of  its  initial  shape;  the 
distance  traveled  is  determined  by  the  group  velocity.  Because  the  disturbance  is 
stationary  in  the  reference  frame  of  the  ship,  the  wake,  which  travels  along  with  it,  is  also 
stationary  in  this  frame.  This  implies  that  each  wave  in  the  wake  must  satisfy  a  phase 
relation  based  on  the  phase  velocity;  the  component  of  the  ship  velocity  along  the  wave 
vector  must  be  equal  to  the  phase  velocity.  Otherwise  a  crest  would  not  be  stationary  in 
this  reference  frame. 

While  the  latter  simplification  involving  the  group  velocity  might  seem  restrictive,  the 
approach  is  related  to  geometrical  optics  and  closely  resembles  the  mathematical  methods 
of  integration  known  as  stationary  phase  and  steepest  descents.  It  is  known  that  these 
methods  frequently  provide  results  that  are  approximately  correct  over  a  domain  that  is 
much  broader  than  might  reasonably  be  expected.  Thus  the  “geometrical”  method 
provides  plausible  results  even  close  to  the  ship. 

The  theory  and  simulation  of  internal  waves  on  a  horizontal  diffuse  layer  can  be 
understood  in  terms  of  a  waveguide  analogy  from  electromagnetic  theory.  In  a  waveguide 
the  waves  within  it  are  bounded  by  the  walls  and  propagate  between  them,  where  they  are 
reflected.  Nevertheless,  the  propagation  can  be  described  by  wave  vectors  parallel  to  the 
waveguide  axis.  At  a  given  frequency  there  may  be  several  possible  wave  vectors  and 
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each  one  corresponds  to  a  “mode”.  Similarly,  at  a  given  frequency,  which  must  be  less 
than  local  B-V  frequency,  waves  can  propagate  with  a  component  in  the  vertical  direction 
but  are  reflected  before  they  reach  a  level  at  which  the  B-V  frequency  equals  the  wave 
frequency.  Therefore  the  waves  tend  to  be  constrained  to  propagate  within  the  layer. 
Outside  these  bounds,  the  wave  amplitudes  fall  off  exponentially  in  a  manner  that 
depends  on  the  wave  frequency,  the  B-V  frequencies  and  how  the  latter  vary  with  depth. 

Wake  simulation  for  the  discrete  internal  layer  has  been  described  in  previous  reports  [2] 
[3].  There  are  significant  differences  between  the  wakes  produced  on  a  discrete  layer 
compared  with  those  on  a  diffuse  layer.  The  discrete  layer  is  equivalent  to  a  very  narrow 
diffuse  layer.  While  the  B-V  frequency  of  a  diffuse  layer  is  usually  limited  to  about  0.1 
rad/s  (as  in  the  Loch  Linnhe  trials)  and  is  usually  much  less,  that  the  B-V  frequency  in  the 
discrete  layer  is  technically  infinite.  These  differences  influence  the  crest  pattern  because, 
as  just  noted,  the  maximum  frequency  of  horizontal  waves  on  a  layer  is  equal  to  the 
maximum  B-V  frequency  of  the  layer.  The  lowest  frequency  waves  in  a  wake  tend  to  lie 
near  the  wake  edge  and  the  highest  frequency  waves  near  the  wake  axis.  Therefore  the 
finite  upper  limit  to  the  B-V  frequency  affects  those  wake  waves  near  the  wake  axis.  The 
effect  on  the  crest  pattern  is  to  shift  the  apex  of  subsequent  crests  backwards,  so  that  there 
is  a  gap  between  each  crest  on  the  axis  of  the  wake.  This  axis  of  course  corresponds  to  the 
ship  track.  The  lower  the  maximum  B-V  frequency,  the  greater  the  backward  shift.  In  the 
case  of  a  discrete  layer  there  is  no  shift  and  all  crests  intersect  at  the  wake  apex. 

Additional  differences  occur  because  waves  can  propagate  with  a  vertical  component. 
This  implies  that  a  disturbance  actually  located  in  the  layer  will  produce  waves  directly 
and  these  can  be  of  large  amplitude  and  can  propagate  vertically  through  the  layer;  this 
tends  to  increase  wake  visibility.  In  comparison  the  disturbance  on  a  discrete  layer  will 
usually  be  located  some  vertical  distance  away  from  it  and  its  effect  at  the  layer  will  be 
attenuated.  Unfortunately  broad  layers  tend  to  exhibit  low  B-V  frequencies  and  the  crest 
pattern  shifts  may  reduce  wake  visibility.  Thus,  as  far  as  detectability  is  concerned,  these 
effects  tend  to  act  in  opposition. 

An  example  of  a  (airborne)  radar  image  of  a  wake  from  the  Loch  Linnhe  trials  is  shown 
in  Figure  1.  The  radar  frequency  is  at  L-band  and  the  ship,  which  was  the  Sir  Tristram 
and  appears  as  the  bright  object  to  the  right,  is  traveling  at  a  speed  of  2  m/s.  The  length  of 
the  image  corresponds  to  a  distance  of  about  4  km  and  the  original  pixel  size  was  about  4 
m. 
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Figure  1.  Sir  Tristram  wake.  Ship  speed  is  2.0  m/s.  From  [5]. 

Internal  Waves 

When  stratification  of  fluid  density  is  present,  internal  waves  can  occur.  These  are  gravity 
waves  and  their  mechanism  resembles  that  of  the  waves  on  the  ocean  surface.  It  is  usually 
assumed  that  the  stratification  is  horizontal  so  that  the  properties  of  the  layer  are 
described  by  vertical  profiles,  such  as  density  and  temperature.  In  fact  the  most  useful 
profile  is  that  of  the  B-V  frequency.  As  noted,  the  B-V  frequency  is  the  frequency  of 
vertical  buoyancy  oscillations  of  a  small  fluid  parcel  that  has  been  disturbed  from  its 
equilibrium  position  and,  in  contrast  to  surface  gravity  waves,  internal  waves  are  three- 
dimensional. 

Oceanic  internal  waves  are  frequently  observed  by  space-bome  Synthetic  Aperture  Radar 
(SAR).  A  useful  description  has  been  provided  by  Jackson  et  al.  (2013)  [7]  as  well  as  in 
[8].  The  waves  are  usually  a  result  of  tidal  flows  over  sills  or  other  obstacles  on  the  ocean 
floor.  They  are  analogous  to  the  waves  on  the  surface  of  a  stream  that  flows  over  a  stone. 
However,  naturally  occurring  internal  waves  that  are  observed  by  SAR  are  typically  non¬ 
linear  wave  packets  and  individual  waves  develop  into  a  soliton.  Internal  wave  solitons 
are  analogous  to  tidal  bores  that  can  occur  on  some  rivers.  In  a  sequence  of  tidal  cycles,  a 
train  of  internal  solitons  can  be  produced.  The  phase  speed  of  the  waves  can  vary 
between  0.3  m/s  to  over  3  m/s;  amplitudes  on  the  layer  can  range  from  a  few  metres  to 
more  than  100  m. 

The  waves  are  visible  in  SAR  imagery  because  surface  flows  are  induced  that  affect  the 
radar  scattering  cross-section  of  the  ocean  [9].  These  flows  compress  or  expand  the  Bragg 
waves  responsible  for  the  radar  backscattering.  In  addition,  periodic  horizontal  flows  can 
create  bands  of  slicks,  which  tend  to  suppress  the  Bragg  waves  [7].  The  waves  can  be 
distinguished  from  other  features  in  radar  imagery  by  the  fact  that  natural  internal  waves 
appear  in  bands  and  the  first  or  leading  bright  or  dark  crest  is  separated  from  the  second 
by  a  slightly  larger  distance  than  that  between  the  second  and  third.  This  is  because  large 
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amplitude  waves  move  faster  than  small  amplitude  waves  and  the  fact  that  the  soliton 
amplitude  tends  to  increase  as  the  wave  packet  evolves  and  leading  waves  coalesce. 

Internal  waves  in  SAR  imagery  tend  to  be  located  on  continental  shelves,  in  straits  and  in 
marginal  seas  [7].  This  is  because  not  only  does  suitable  horizontal  stratification  but  also 
the  conditions  for  wave  generation  exist.  It  does  not  imply  that  layers  capable  of 
supporting  internal  waves  are  not  present  in  the  open  ocean.  However,  density  variations 
associated  with  fresh  water  flowing  over  salt  water  will  tend  to  be  confined  to  coastal 
regions  and  the  arctic  (ice  melt).  In  the  open  ocean,  temperature  variations  are  expected 
to  be  the  main  factor.  Nevertheless,  the  global  distribution  of  natural  internal  waves  may 
be  useful  for  estimating  the  likelihood  of  a  significant  layer  in  a  given  area  of  ocean. 

Global  Ocean  Associates  has  created  an  atlas  of  internal  wave  occurrences  [10].  This 
contains  a  chapter  by  Apel  summarizing  the  background  and  theory  of  the  waves.  He 
remarks  (in  2004)  that  the  distribution  of  deep  water  internal  waves  is  still  a  research 
project  and  this  seems  to  be  confirmed  in  [7],  The  atlas  contains  details  about  the 
generation  and  interpretation  of  internal  waves  and  the  section  on  the  Strait  of  Gibraltar  is 
particularly  interesting.  The  flow  of  the  upper  layer  of  ocean  water  is  eastwards  from  the 
Atlantic  Ocean  through  the  Strait.  The  lower  layer  of  denser  water  flows  westward  from 
the  Mediterranean.  Internal  waves  are  generated  by  these  flows  over  very  pronounced 
bottom  features,  such  as  the  Camarinal  Sill.  Waves  are  present  year-round  particularly  on 
the  eastern  side  of  the  Strait.  In  summer,  waves  can  be  produced  by  westward  flows  on 
the  western  side  due  to  the  presence  of  a  strong  thermocline  in  the  Atlantic.  The  waves 
are  strongest  during  the  spring  tide  and  weakest  at  neap  tide.  This  effect  is  associated 
with  the  strength  of  the  flows  through  the  Strait. 

An  example  of  a  B-V  profile  for  the  Strait  of  Gibraltar  is  provided  in  [10]  and  this  is 
adopted  for  one  canonical  model  in  this  report.  Figure  2,  from  [8],  shows  the  distribution 
of  internal  wave  observations  over  a  period  from  2002  to  2004. 
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Figure  2.  Locations  of  3500  occurrences  of  non-linear  internal  waves 
In  MODIS  true-color  images  from  2002  to  2004  (from  [8]). 
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In  summary,  significant  internal  layers  can  certainly  be  found  on  most  if  not  all 
continental  shelves,  especially  in  summer.  They  can  also  be  found  in  the  open  ocean 
especially  in  the  tropics  during  the  summer. 

Layer  Models 

Two  basic  models  are  employed  for  the  vertical  B-V  frequency  profile.  The  first  relies  on 
a  Lorentzian  distribution  with  depth.  The  parameters  are  the  maximum  B-V  frequency, 
the  depth  of  this  peak  and  its  width  and  the  depth  of  the  ocean.  An  example  is  shown  in 
Fig.  3,  which  represents  the  profile  off  Gibraltar.  The  peak  B-V  frequency  is  0.015  rad/s 
at  a  depth  of  60  m.  The  width  is  100  m.  There  is  no  mixed  layer  at  the  surface. 

The  second  type  of  profile  is  read  from  a  file  and  is  believed  to  be  representative  of  the 
profile  during  the  Loch  Linnhe  trials.  This  is  derived  from  Watson  et  al.  [5]  and  is  shown 
in  Fig.  4.  The  layer  is  close  to  the  surface  and  is  quite  narrow;  the  maximum  B-V 
frequency  is  correspondingly  greater  than  that  of  the  Gibraltar  profile.  Again  there  is  no 
significant  mixed  layer  near  the  surface.  As  will  be  shown,  a  thick  mixed  layer  tends  to 
reduce  the  surface  velocities  and  the  likelihood  of  observing  the  wake. 
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Figure  3.  Gibraltar  vertical  profile. 
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Figure  4.  Loch  Linnhe  vertical  profile. 


Several  other  models  are  available  from  [10].  These  can  be  represented  approximately  by 
a  Lorentz  distribution  with  the  parameters  in  Table  1. 


Table  1. 

Internal  Layer  Parameters 


Area 

Peak 

N 

(rad/s) 

Peak 

Depth 

(m) 

Width 

(m) 

Water 

Depth 

(m) 

Mixed 

Depth 

(m) 

Celtic  Sea 

0.015 

30 

50 

200 

0 

Gulf  of  Mexico 

0.045 

15 

20 

275 

5 

New  England  Shelf 

0.052 

15 

10 

75 

0 

New  York  Bight 

0.070 

15 

10 

300 

0 

South  West  Africa 

0.025 

25 

30 

250 

0 

Sulu  Sea 

0.020 

120 

200 

1800 

0 

US  Pacific  (N.W.) 

0.090 

6 

6 

30 

2 

Andaman  Sea 

0.023 

100 

150 

1000 

0 

Hull  Models 

A  number  of  hull  models  are  available  for  the  simulation.  The  primary  models  are  the 
Taylor  Standard  Series,  which  are  based  on  a  British  cruiser  and  are  applicable  to  twin 
screwed  warships  [11],  and  the  David  Taylor  Model  Basin  (DTMB)  Series  60  models 
that  are  applicable  to  single  screwed  merchant  ships  [12].  In  [3]  it  was  demonstrated  that, 
other  than  the  principal  dimensions,  the  most  important  hull  parameter  is  the  block 
coefficient,  C«,  which  represents  the  ratio  of  the  volume  of  the  submerged  hull  to  the 
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volume  of  a  rectangular  parallelpiped  that  just  encloses  the  submerged  hull.  While  the 
Taylor  model  can  be  implemented  over  a  continuous  range  of  block  coefficients,  the 
DTMB  model  includes  5  separate  hull  designs  covering  a  range  0.6  <  Cg  <  0.8.  Warships 
designed  for  speed  are  generally  slender  with  a  low  coefficient  but  merchantmen 
generally  he  at  the  upper  end  of  the  range. 

Secondary  models  are  the  Wigley  model  [13]  for  which  Cg  =  4/9,  a  single  source 
separated  from  a  single  sink  (doublet)  model  and  a  single  source. 

The  ship  models  are  published  as  tables  of  hull  offsets  or  transverse  coordinates  at 
various  positions  or  “stations”  along  the  length  of  the  hull.  The  data  points  tend  to  be 
widely  spaced  and  intermediate  offsets  and  longitudinal  positions  must  be  interpolated; 
numerical  splines  are  the  preferred  choice  as  this  mimics  the  procedures  that  were 
formerly  used  in  ship-building.  The  Taylor  series  is  based  on  a  parent  form  (a  plan  view 
from  keel  to  waterline  is  shown  in  Figure  5)  and  an  algorithm  developed  by  Gertler 
allows  its  generalization  into  a  series.  Unfortunately  some  tables  contain  obvious 
typographical  errors  and  in  the  Taylor  Series  there  are  some  other  minor  errors  needing 
correction  [14].  Also  the  algorithm  is  only  satisfactory  for  0.5  <  Cg  <  0.7.  Outside  this 
range,  the  hull  forms  become  increasingly  unrealistic. 

A  similar  plot  of  a  DTMB  form  using  splines  is  shown  in  Figure  6.  These  plots 
demonstrate  that  spline  interpolation  is  satisfactory. 

Figures  7  to  9  show  plots  of  the  offsets.  The  bow  section  offsets  are  on  the  right  hand  side 
of  the  plot  and  the  stem  section  offsets  are  on  the  left.  “FP”  and  “AP”  are  the  forward  and 
aft  perpendiculars  respectively.  “WL”  stands  for  waterline.  The  waterline  for  the  Taylor 
model  corresponds  to  the  horizontal  line  at  about  half  the  height  of  the  plot.  For  the 
DTMB  model,  it  lies  at  “1.00WL”  and  for  the  Wigley  model  it  corresponds  to  the  lower 
horizontal  line.  The  numbers  are  labels  for  the  longitudinal  stations.  It  is  worth  noting  the 
small  bulbous  bow  on  the  Taylor  series  and  its  absence  on  the  DTMB  series.  The 
submerged  aft  section  on  the  Taylor  series  is  narrower  than  on  the  DTMB  series  to 
accommodate  and  support  twin  screws.  The  Wigley  hull  is  parabolic  in  the  transverse  and 
longitudinal  directions  but  the  draft  is  constant  unlike  a  true  canoe  shape.  The  example  in 
Figure  9  also  has  freeboard. 
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Distance  from  Bow 

Figure  5.  Plan  view  of  Taylor  parent  hull  offsets  over  150  stations  derived  from  original 
29  stations  using  a  cubic  spline  interpolation;  C«  ~  0.493.  Red  to  light  blue  lines  represent 
offsets  at  equal  vertical  intervals  from  close  to  keel  to  waterline. 


Figure  6.  Plan  view  of  DTMB  Hull  offsets  over  150  stations  derived  from  original  25 
stations  using  a  cubic  spline  interpolation;  Cb  =  0.7.  Blue  to  purple  lines  represent  offsets 

at  equal  intervals  from  keel  to  waterline. 


Figure  7.  Example  of  Taylor  series  hull  offsets  (from  [1 1]). 
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Figure  9.  Example  of  Wigley  hull  offsets  (from  [13]). 


Figure  10.  Source  distribution  for  Wigley  hull  over  150  stations.  The  lines,  blue  to  purple 
represent  the  distributions  from  the  keel  to  the  waterline  at  equal  intervals  of  the  vertical 
coordinate.  The  hull  beam  and  length  are  one  unit  of  length  and  the  source  strength  is  a 

volume  per  unit  time. 
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Figure  11.  Source  distribution  for  Taylor  hull  (Cs  =  0.7)  over  150  stations.  The  lines,  blue 
to  purple  represent  the  distributions  from  the  keel  to  the  waterline  at  equal  intervals  of  the 
vertical  coordinate.  The  hull  beam  and  length  are  one  unit  of  length  and  the  source 

strength  is  a  volume  per  unit  time. 


11 


Figure  12.  Source  distribution  for  DTMB  Series  60  hull  (Cg  =  0.7)  over  150  stations.  The 
lines,  blue-red-green-purple,  represent  the  distributions  from  the  keel  to  the  waterline  at 
equal  intervals  of  the  vertical  coordinate.  The  hull  beam  and  length  are  one  unit  of  length 
and  the  source  strength  is  a  volume  per  unit  time. 

Figures  10  to  12  show  the  source  distribution  curves  for  the  Wigley,  Taylor  and  DTMB 
Series  60  hulls.  In  the  last  two  cases  the  block  coefficient  is  0.7.  The  distributions  over  4 
horizontal  slices  are  shown  at  equal  height  intervals.  While  the  distribution  for  the 
Wigley  hull  is  quite  simple,  the  distributions  for  the  practical  hulls  are  much  more 
complicated.  The  bow  and  stem  sections  are  separated  and  there  is  a  distinct  midships 
section  where  the  distribution  is  close  to  zero.  The  distribution  of  the  bow  section  reflects 
efforts  to  improve  seakeeping  and  reduce  deck-wetting  as  well  as  to  reduce  wave 
resistance;  the  stern  distribution  reflects  the  necessity  of  accommodating  either  twin 
screws  and  shafts  or  a  single  screw. 

Theory 

This  is  based  on  previous  work  [1]  and  a  summary  is  as  follows. 


In  principle,  waves  on  the  layer  at  a  given  depth  are  represented  by  a  sum  of  horizontally 
propagating  waves.  The  amplitudes  of  a  wave  as  a  function  of  depth  are  solutions  of  a 
differential  equation  subject  to  boundary  conditions  that  must  be  solved  numerically: 


d2Q-_ 

dz 2 


+  k~ 


N\z) 


-1 


CO 


Q-=  0 


(1) 


Here  Qz  is  the  Fourier  transform  of  the  vertical  flux,  N  is  the  B-V  frequency,  co  is  the 
angular  frequency  and  k  is  the  angular  wave  number.  The  equation  is  sometimes  known 
as  the  Taylor-Goldstein  equation  but  in  our  case  it  is  the  simplest  version  because 
currents  and  current  shears  are  neglected.  The  boundary  conditions  are  simplified  by 
neglecting  the  usual  surface  gravity  waves  but  this  has  no  practical  effect  on  the  solution 
because  the  vertical  displacements  due  to  internal  waves  are  very  small.  The  vertical  fluid 
velocity  is  set  to  zero  at  the  fluid  surface  as  well  as  at  its  floor.  Because  wake  waves  exist 
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over  a  wide  range  of  frequencies,  the  differential  equation  must  be  solved  over  this  range. 
At  each  frequency,  a  solution  of  the  equation  yields  a  sequence  of  eigenvalues  and 
eigenvectors.  Each  member  of  the  sequence  represents  a  mode  of  oscillation.  An 
eigenvalue  is  the  angular  wave  number  of  the  horizontal  wave  and  the  corresponding 
eigenvector  is  its  nonnalized  amplitude  as  a  function  of  depth.  In  principle,  the  final 
amplitude  is  a  sum  over  these  eigenvectors;  the  coefficients  would  depend  on  how  the 
waves  are  generated. 

The  Taylor-Goldstein  equation  is  of  the  Sturm-Liouville  type.  However,  the  weighting 
factor,  conventionally  denoted  by  p,  is  not  always  positive  as  is  usually  assumed  in 
standard  texts,  such  as  [15]  or  [16].  (The  situation  is  further  confused  by  the  omission  of 
the  factor  p  in  an  equation  at  the  top  of  page  295  in  [15].)  This  is  because  the  application 
is  usually  to  waves  on  non-homogeneous  strings  where  the  solutions  are  always 
oscillatory.  In  the  present  application  to  internal  waves,  the  waves  are  constrained  within 
the  layer  and  tend  to  fall  off  exponentially  above  and  below  it;  in  these  regions  p  is 
negative.  Therefore  it  seems  prudent  to  confirm  that  the  usual  theory  is  applicable.  The 
theory  is  discussed  in  Appendix  1 ,  where  it  is  shown  that  the  only  potential  problem  is  a 
failure  to  nonnalize  the  eigenvector  at  some  frequencies.  It  is  worth  noting  that  so  far  this 
has  never  occurred  numerically  and  it  may  reasonably  be  concluded  that  the  usual  Sturm- 
Liouville  theory  is  satisfactory. 

The  lowest  mode  of  oscillation  is  a  sinuous  mode  in  which  the  entire  layer  moves  up  and 
down  with  amplitudes  that  depend  on  depth  and  are  proportional  to  a  sum  of  the 
eigenvectors  of  this  lowest  mode.  Clearly  these  eigenvectors  have  no  zero-crossings.  If  a 
vessel  is  above  or  below  the  layer  peak,  as  it  moves  horizontally,  it  tends  to  push  the  layer 
down  or  up  and  the  lowest  mode  will  predominate. 

When  a  vessel  is  close  to  the  peak  of  the  layer,  it  tends  to  push  the  layer  more  or  less 
equally  in  both  directions.  This  corresponds  to  a  varicose  mode  in  which  the  propagating 
waves  cause  the  layer  to  vary  in  thickness.  This  corresponds  to  the  next  lowest  mode  and 
its  eigenvectors  have  a  single  zero-crossing.  It  turns  out  that  the  propagation  velocities  of 
the  waves  decrease  rapidly  as  the  mode  number  increases  so  that  the  lowest  mode 
contributes  the  broadest  wake  and  is  the  most  important  by  far;  higher  modes  can  usually 
be  neglected.  The  next  mode  is  important  only  when  a  vessel  is  traveling  close  to  the 
layer  peak.  This  situation  is  not  considered  in  detail  here. 

The  relation  between  the  frequency  and  the  horizontal  wave  number  in  a  given  mode 
represents  a  dispersion  relation.  This  is  conveniently  calculated  by  solving  the  Taylor- 
Goldstein  equation  at  equal  intervals  of  the  wave  frequency  up  to  the  maximum  B-V 
frequency.  The  wave  number  at  intermediate  frequencies  can  be  found  by  linear 
interpolation.  As  described  in  Appendix  2,  the  dispersion  equation  can  be  employed  to 
find  the  crest  pattern  in  the  vessel  wake.  The  crest  pattern  for  any  mode  does  not  depend 
on  the  depth  of  the  source  generating  the  waves  though  this  does  affect  the  amplitude  of 
waves  in  a  wake  [1],  The  dispersion  relation  for  the  Gibraltar  profile  is  shown  in  Fig.  13. 
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Figure  13.  Dispersion  relation  of  the  lowest  mode  for  the  Gibraltar  profile. 

At  very  low  wave  numbers,  the  graph  is  linear  and  its  slope  is  equal  to  both  the  group  and 
phase  velocities.  As  expected,  the  maximum  frequency  is  equal  to  the  maximum  B-V 
frequency.  Dispersion  relations  for  all  reasonable  profiles  have  similar  characteristics  and 
resemble  this.  The  dispersion  relation  for  the  Loch  Linnhe  profile  is  provided  in  Fig.  14. 
The  graph  exhibits  some  additional  minor  structure  (at  about  2.0  rad/m)  associated  with 
the  complex  structure  in  the  original  profile. 


Figure  14.  Dispersion  relation  of  the  lowest  mode  for  the  Loch  Linnhe  profile. 

The  eigenfunctions,  which  depend  on  frequency,  represent  the  variation  of  amplitude  of  a 
mode  at  a  given  frequency.  At  low  frequencies  they  are  easily  found  during  the  solution 
of  the  Taylor-Goldstein  equation.  At  high  frequencies,  close  to  the  maximum  B-V 
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frequency,  numerical  problems  can  occur  and  these  are  discussed  later.  Examples  are 
shown  for  the  first  two  modes  at  intermediate  frequencies  in  each  of  the  two  layer  profiles 
in  Figs.  15  to  18. 


Figure  15.  Eigenfunction  for  lowest  mode  of  Gibraltar  profile  for  frequency  of  0.0075 
rad/s.  The  amplitude,  which  is  in  arbitrary  units,  is  proportional  to  the  vertical  flux. 
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Figure  16.  Eigenfunction  for  next  lowest  mode  of  Gibraltar  profile  for  frequency  of 

0.0075  rad/s. 
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Figure  17.  Eigenfunction  for  lowest  mode  of  Loch  Linnhe  profile  for  frequency  of  0.0625 

rad/s. 


Figure  18.  Eigenfunction  for  next  lowest  mode  of  Loch  Linnhe  profile  for  frequency  of 

0.0625  rad/s. 

In  the  lowest  mode  the  low-frequency  eigenfunctions  seem  to  be  a  smoothed  version  of 
the  B-V  profile.  At  higher  frequencies,  co,  close  to  the  maximum  B-V  frequency, 
eigenfunctions  peak  strongly  when  the  B-V  frequency,  N  >  co.  The  Loch  Linnhe  profile 
exhibits  two  peaks  here  and  indeed  there  can  be  two  smoothed  peaks  in  the  eigenfunction. 

Waves  in  a  wake  are  produced  as  the  vessel  creates  a  disturbance  in  the  water  during  its 
movement,  which  for  our  purposes  is  a  constant  horizontal  velocity.  Waves  propagate 
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outwards  from  their  point  of  creation  at  the  group  velocity.  At  a  given  point  in  the  wake, 
which  is  stationary  in  the  reference  frame  of  the  vessel,  waves  interfere  constructively  so 
that  a  phase  condition  must  be  met.  This  is  equivalent  to  insisting  that  the  component  of 
the  vessel  velocity  in  the  direction  of  the  wave  must  be  equal  to  the  phase  velocity  of  the 
wave.  If  the  phase  of  the  wave  is  defined  to  be  zero  at  the  ship,  the  first  crest  occurs  when 
the  phase  is  2n,  etc.  Using  the  dispersion  curve,  such  as  that  in  Figure  14,  these  concepts 
can  be  combined  and  are  sufficient  to  calculate  the  crest  pattern.  The  method  is  described 
in  Appendix  2,  where  it  is  shown  that  the  locus  of  a  crest  is  determined  in  terms  of  the 
phase  and  group  velocities  and  how  they  vary  with  frequency. 

For  example,  the  pattern  for  the  first  5  crests  is  shown  in  Figure  19  for  a  ship  traveling  at 
2  m/s  over  the  profile  in  Figure  4.  This  can  be  compared  to  the  image  of  a  wake  from  the 
Loch  Linnhe  trials  in  Figure  1.  The  pattern  is  a  good  match  to  the  radar  observations.  The 
crest  pattern  for  the  next  mode  is  similar  except  that  the  wake  is  much  narrower. 


Figure  19.  The  lowest  mode  crest  pattern  for  the  Loch  Linnhe  profile  showing  the  first  5 
crests,  which  are  the  most  prominent.  The  ship  speed  is  2  m/s  and  the  image  length  is  4 
km.  The  ship  is  located  at  the  right  hand  edge  of  the  picture  box. 
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In  [1]  the  theory  was  sketched  and  (inserting  the  correct  coefficient  for  the  inverse 
transforms)  it  can  be  shown  that  the  vertical  component  of  the  flux,  qz,  for  a  stationary 
oscillating  source  at  depth  c  is  given  by: 

1  fdF,  ' 

(2tt)3  J  4"  k2-k2 


q,  = 


dz 


i(«- k.r) 


dkdco 


(2) 


2z=c 


Here,  z  is  the  vertical  coordinate  (positive  downwards),  k„  and  Q„  are  the  eigenvalue  and 
the  orthonormal  eigenvector  for  mode  n;  co  and  k  are  the  wave  frequency  and  wave 
vector,  respectively.  The  source  strength  is  So  and  this  can  be  expressed  as  a  mass  of  fluid 
per  unit  time.  The  function  F„  is  given  by: 


F,=QnV) 
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where  N  is  the  B-V  frequency.  If  we  are  only  interested  in  the  lowest  mode,  which  we 
label  ‘O’,  and  perform  the  integral  over  frequency,  we  have: 
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We  now  shift  the  frame  of  reference  from  the  fluid  to  a  constant  source  moving  at 
velocity  U  and  combine  this  with  the  Doppler  equation,  which  is: 


0  =  0)n  -  k  ■  U 


(5) 


Then  we  have  for  the  wake: 
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(6) 


where 

r'  =  r-Ut  (7) 

Clearly  r'  is  the  position  of  a  point  in  the  wake  relative  to  the  source.  A  dimensional 
analysis  can  be  implemented  using  Table  2  (this  includes  variables  used  in  [1]);  it  is 
easily  verified  that  (6)  is  dimensionally  correct. 


Table  2. 

Variable  Dimensions 


Variable 

Dimensions 

q 

M/(L2T) 

k 

1/L 

dk 

1/L2 

CO 

1/T 

s 

M/(LJT) 

So 

M/T 

Q 

M 

On 

~m7n 

As  shown  in  Appendix  2,  at  each  point  in  the  wake  a  wave  of  specific  frequency  and 
wave  number  predominates.  The  values  of  these  parameters  can  be  established 
geometrically.  This  would  allow  the  integral  in  (6)  to  be  evaluated  by  determining  the 
value  of  ko  at  a  point.  The  phase  of  the  wave  is  also  detennined. 
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Appendix  3  provides  the  details  of  the  detennination  of  the  horizontal  surface  flow  and 
the  evaluation  of  the  integral  using  the  stationary  phase  method.  The  horizontal  flow  is 
related  to  qz  by  the  fact  that  the  divergence  of  the  velocity  is  zero  (water  is 
incompressible).  The  limits  of  integration  for  k  can  be  extended  over  [-00,  00]  because 
waves  in  the  far  wake  must  propagate  at  an  angle  9  to  the  ship’s  track  and  cannot 
propagate  in  the  reverse  direction.  They  can  also  be  considered  to  be  slightly  damped  to 
establish  that  the  pole  that  contributes  to  the  integral  in  a  contour  integration  lies  in  the 
upper  half  of  the  complex  plane.  The  remaining  integral  over  9  is  evaluated  by  a 
stationary  phase  technique  as  was  employed  by  Tunaley  [17].  Subsequently  Chunchuzov 
et  al.  [18]  appear  to  have  used  a  similar  approach. 

Appendix  5  provides  some  comparisons  between  the  discrete  and  diffuse  theories  and 
confirms  consistencies  between  the  two  approaches. 

Numerical  Analysis 

The  first  part  of  the  simulation  is  the  calculation  of  the  eigenvalues  and  eigenvectors  of 
the  Taylor-Goldstein  equation.  The  eigenvalue  represents  the  wave  number  of 
horizontally  propagating  waves  and  its  eigenvector  represents  the  amplitude  of  the 
vertical  component  of  the  mass  flux  as  a  function  of  depth.  This  calculation  must  be 
completed  for  a  range  of  frequencies. 

The  amplitude  of  the  vertical  component  of  the  mass  flux  (fluid  density  times  the 
velocity)  is  zero  at  the  surface  but  its  derivative  with  depth  here  can  be  given  a  small 
positive  value.  If  some  arbitrary  value  is  assigned  to  the  horizontal  wave  number,  the 
differential  equation  can  then  be  marched  forward.  Typically,  the  amplitude  will  increase 
with  depth  as  the  maximum  B-V  frequency  is  reached  and  then  it  will  either  increase  or 
decrease  rapidly  reaching  large  positive  or  negative  values  at  the  bottom  of  the  fluid.  The 
desired  solution  is  the  value  of  the  wave  number  when  the  amplitude  at  the  bottom  is 
zero. 

One  way  to  achieve  this  is  to  perform  the  integration  using  the  Runge-Kutta  (R-K) 
algorithm.  Implicit  in  this  is  the  idea  that  the  B-V  frequency  is  reasonably  smooth.  This 
algorithm  appears  to  function  successfully  at  frequencies  up  to  about  %  of  the  maximum 
B-V  frequency  and  delivers  reasonable  eigenvalues  and  eigenvectors  but  then  several 
problems  start  to  emerge.  These  involve  the  precision  of  the  floating  point  numbers  and 
the  accuracy  of  the  R-K  algorithm. 

The  eigenvalue  can  be  established  using  a  bisection  search  to  locate  the  appropriate 
horizontal  wave  number  to  make  the  amplitude  at  the  bottom  zero.  As  the  frequency 
increases  beyond  a  critical  fraction  of  the  B-V  frequency,  the  sensitivity  of  this  amplitude 
grows  to  such  an  extent  that  double  precision  is  wholly  insufficient  to  locate  the 
eigenvalue  such  that  the  eigenvector  is  satisfactory.  Therefore  an  eigenvalue  can  be 
obtained  accurately  but  the  eigenvector  may  exhibit  a  sharp  upturn  or  downturn  of  large 
amplitude  near  the  ocean  bottom.  In  extreme  cases,  the  eigenvector  is  useless.  Even 
though  the  code  is  in  C-Sharp,  which  is  very  forgiving  of  floating  point  overflow,  “Not  a 
Number”  (NaN)  was  returned  for  some  calculations.  In  part,  these  effects  may  be 
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attributed  to  small  errors  introduced  by  the  R-K  algorithm  itself.  Therefore  the  R-K 
algorithm  was  replaced  by  a  matrix  approach,  which  is  described  in  Appendix  4. 

Unlike  the  R-K  model,  the  model  of  the  layer  in  the  matrix  approach  is  stepwise  in 
nature,  so  that  the  B-V  frequency  can  change  discretely  from  step  to  step  but  is  constant 
over  each  step.  The  solution  of  the  Taylor-Goldstein  equation  is  either  sinusoidal  (if  the 
frequency  is  less  than  the  local  BV  frequency)  or  exponential  (if  it  is  greater).  Any 
second  order  differential  equation  can  be  replaced  by  two  single  order  differential 
equations  so  that  a  marching  solution  can  be  expressed  as  a  matrix  multiplication  with  a 
two-component  vector.  The  vector  components  are  typically  the  amplitude  and  its 
derivative  (see  Appendix  4).  Therefore,  subject  only  to  rounding  errors  and  the 
limitations  of  the  stepwise  model,  the  integration  can  be  performed  exactly. 

When  this  model  was  implemented,  the  problems  of  the  limited  precision  and  the 
sensitivity  of  the  eigenvector  to  the  eigenvalue  remained.  However,  by  limiting  the 
magnitude  of  the  variables,  the  NaN  returns  can  be  eliminated  and  the  upturns  or 
downturns  can  simply  be  set  to  zero  (or  just  ignored).  This  provides  satisfactory 
eigenvectors  over  the  entire  range  of  frequencies  and  wave  numbers  that  are  required. 

The  present  algorithm  functions  well  for  the  Loch  Linnhe  and  Gibraltar  profiles  and 
yields  close  but  not  exactly  similar  eigenvalues  and  eigenvectors  to  those  using  the  R-K 
approach.  Small  discrepancies,  such  as  those  observed,  are  expected  because  the  models 
are  not  identical. 

While  the  development  was  principally  for  the  lowest  mode,  the  algorithm  was  tested  for 
the  second  mode,  which  has  a  single  zero-crossing.  It  functions  well  for  this  mode  too. 

The  Program 

The  program  was  constructed  using  Microsoft’s  Visual  Studio  in  the  C-Sharp  language. 
The  main  window  is  shown  in  Figure  20.  The  B-V  profile  can  be  loaded  from  a  file  or  be 
created  using  a  Lorentzian  model  by  clicking  on  the  appropriate  radio  button.  In  the  latter 
case  the  B-V  frequency  is  finite  at  the  surface  and  there  is  no  mixed  layer,  which  is 
similar  to  the  standard  Loch  Linnhe  profile. 

The  operator  must  also  specify  the  length  of  the  wake  to  display  and  an  offset  for  the 
velocity  profile  along  a  track  that  is  parallel  to  the  wake  axis. 

There  is  a  choice  of  hull  type  using  radio  buttons  and  the  vessel’s  principal  dimensions 
can  be  entered.  The  “Go”  button  is  then  clicked.  If  the  hull  model  is  a  Taylor  or  DTMB 
type,  a  subsidiary  window  allows  the  operator  to  choose  a  block  coefficient  within  limits 
or  from  a  list,  respectively.  When  the  hull  is  a  single  source,  a  doublet  or  submarine,  a 
subsidiary  window  allows  the  operator  to  enter  a  depth.  The  source  strengths  for  these  are 
calculated  from  the  vessel’s  draft,  T,  and  beam,  B ,  using  an  elliptical  transverse  profile 
model.  Thus  the  source  strength,  So  =  nBTU/4  (U  is  the  source  speed).  The  separation  of 
the  doublet  (a  source  of  strength  So  and  an  equal  and  opposite  sink)  is  set  to  the  vessel’s 
length,  L. 
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The  output  is  a  display  of  the  surface  velocity  component  transverse  to  the  axis  of  the 
wake.  Beneath  this  is  a  graph  of  the  transverse  velocity  component  along  a  track  that  is 
parallel  to  the  wake  axis  and  offset  from  it  by  the  specified  distance.  The  time  on  a 
typical  personal  computer  running  Microsoft’s  Windows  operating  system  is  usually  less 
than  20  s. 


□9  Diffuse  Internal  Layer  Simulation 
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Figure  20.  Main  window  of  simulation  program. 

The  example  in  Figure  20  is  for  the  Loch  Linnhe  B-V  profile  and  the  Sir  Tristram  ship. 
The  hull  type  is  Wigley.  Differences  in  the  appearance  of  the  wake  in  the  main  display  do 
occur  when  the  ship  hull  is  changed  to  a  Taylor  or  DTMB  model  and  when  the  block 
coefficient  is  changed.  However,  these  are  quite  minor  and  are  not  easily  discerned. 
However,  they  are  usually  quite  obvious  in  the  graph  beneath. 

The  full  wake  angle  and  the  maximum  transverse  velocity  along  the  offset  track  are  also 
output. 

Simulation  Results 

To  verify  the  current  simulation  of  a  diffuse  layer  and  also  the  previous  simulation  using 
a  discrete  layer,  simulations  were  compared  to  each  other  and  to  the  Loch  Linnhe  trials 
data.  Figure  21  shows  an  example  of  the  transverse  wake  surface  velocity  component  that 
was  derived  from  observational  data  for  a  wake  generated  by  the  Sir  Tristram  landing 
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ship.  These  data  were  band-passed  fdtered  to  reduce  the  noise  components.  The 
effectiveness  of  this  and  the  magnitude  of  the  errors  are  not  clear. 


Figure  21.  Measured  cross-track  current  (and  strain  rate)  corresponding  to  Sir  Tristram 

simulation,  U=  2.0  m/s  (from  [5]). 

Figure  22  shows  a  simulated  version  assuming  the  depth  of  a  discrete  layer  was  3  m.  The 
number  of  oscillations  is  too  large  compared  with  the  observations  but  their  amplitudes 
appear  to  have  the  appropriate  magnitude.  If  the  depth  of  the  layer  is  increased  to  5  m,  the 
number  of  oscillations  is  reduced  but  their  amplitudes  are  too  small. 
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Figure  22.  Simulated  cross-track  surface  velocity  component  using  a  Wigley  hull  and 
depth  of  discrete  interface  of  3  m  (from  [3]). 
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The  corresponding  velocity  profile  for  the  Sir  Tristram  simulation  is  shown  in  Figure  23, 
which  is  similar  in  part  to  the  graph  in  Figure  20.  The  B-V  profile  in  Figure  4  was 
employed  and  it  was  assumed  that  the  Sir  Tristram  had  a  length  of  138  m,  a  beam  of  17  m 
and  a  draft  of  3.9  m  [3];  its  speed  was  2  m/s.  The  resolution  in  depth  for  the  simulation 
was  0.1  m  and  the  horizontal  resolution  was  about  6  m.  The  horizontal  axis  is  in  distance 
rather  than  time  (as  in  Figure  21);  the  conversion  involves  just  a  factor  of  2.  The  number 
of  oscillations  is  of  the  correct  order  and  apparently  the  amplitude  is  greater  than  that 
observed  by  about  30%.  In  addition,  the  general  shape  of  the  profile  is  similar  to  that  of 
the  discrete  simulation  including  the  irregularity  or  shoulder  at  about  1200  m.  The 
agreement  suggests  that  the  simulation  is  providing  reasonable  results  and  that,  unlike  the 
discrete  layer  model,  no  corrections  are  needed  to  achieve  approximate  agreement 
simultaneously  between  the  amplitudes  and  the  number  of  oscillations. 
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Figure  23.  Simulation  of  velocity  profile  169  m  from  ship  track  using  diffuse  layer  and 

Wigley  hull. 

The  effect  of  changing  the  hull  design  to  a  Taylor  or  a  DTMB  Series  60  type  is  shown  in 
Figures  24  to  27. 
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Figure  24.  Simulated  velocity  profile  using  the  previous  parameters  with  a  Taylor  hull 

with  block  coefficient,  Cg,  of  0.5. 
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Figure  25.  Simulated  velocity  profile  using  the  previous  parameters  with  a  DTMB  Series 

60  hull  with  block  coefficient,  Cg,  of  0.6. 


0.08 
0.06 
^  0.04 
0.02 
2“  0 
O  -0.02 

£!  -0.04 
-0.06 
-0.08 

0  1000  2000  3000  4000  5000  6000 


Distance  (m) 


Figure  26.  Simulated  velocity  profile  using  the  previous  parameters  with  a  Taylor  hull 

with  block  coefficient,  Cg,  of  0.  7. 
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Figure  27.  Simulated  velocity  profile  using  the  previous  parameters  with  a  DTMB  Series 

60  hull  with  block  coefficient,  Cg,  of  0.8. 


As  with  the  discrete  layer,  the  main  differences  associated  with  the  hull  are  in  the  leading 
cycles  of  the  velocity  profiles.  These  change  shape  and  the  shoulder  tends  to  become 
more  pronounced  as  the  block  coefficient  rises.  Otherwise  the  amplitudes  of  the  surface 
velocities  remain  much  the  same. 


A  plot  of  the  maximum  transverse  surface  flow  velocity  along  the  parallel  track  against 
ship  speed  is  shown  in  Figure  28  for  the  Sir  Tristram.  The  layer  profile  is  appropriate  to 
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the  Loch  Linnhe  trials  and  the  Taylor  model  was  used  with  a  block  coefficient  of  0.7. 
When  the  ship  speed  is  greater  than  about  2.5  m/s,  the  maximum  transverse  flow  velocity 
barely  changes.  As  the  ship  speed  increases,  the  decrease  in  flow  velocity  is  due  to  the 
additional  distance  of  the  sensor  from  the  ship  as  it  encounters  the  wake.  The  general 
result  seems  to  be  in  accord  with  remarks  in  [5],  in  which  it  is  noted  that  the  maximum 
observed  flow  velocities  from  the  Sir  Tristram  are  always  about  3  cm/s  . 
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Figure  28.  Variation  of  maximum  surface  velocity  with  ship  speed. 


A  mixed  layer  has  a  profound  effect  on  the  amplitudes  of  the  internal  wave  wake.  The 
main  widow  for  a  simulation  with  a  mixed  layer  4  m  deep  above  the  standard  Loch 
Linnhe  profile  is  shown  in  Figure  29;  the  hull  is  of  the  Wigley  type.  The  B-V  frequency 
is  zero  in  this  mixed  layer.  The  Sir  Tristram  draft  is  similar  to  the  layer  depth  so  that  the 
ship  does  not  penetrate  the  layer  where  the  B-V  frequency  is  finite.  There  are  several 
obvious  effects.  Firstly,  the  maximum  amplitude  of  the  surface  velocity  along  the  parallel 
track  is  much  smaller  than  in  the  absence  of  a  mixed  layer.  Secondly,  the  polarity  of  the 
surface  velocity  is  reversed  because  all  sources  over  the  ship  are  now  pushing  the  layer 
downwards  rather  than  the  previous  net  effect  of  an  upwards  push.  Thirdly,  the  wake 
angle  has  increased. 
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Figure  29.  Main  window  of  the  simulation  using  the  Sir  Tristram  ship  with  a  4  m  mixed 
layer  above  the  standard  Loch  Linnhe  profile. 

The  velocity  profile  for  the  above  simulation  is  plotted  in  detail  in  Figure  30.  This  profile 
is  more  complicated  than  those  considered  previously  mainly  because  there  is  now 
significant  interference  between  sources  at  the  bow  and  sinks  at  the  stern.  The  maximum 
transverse  surface  velocity  along  the  parallel  track  is  about  5  mm/s  compared  with 
velocities  of  30  mm/s  during  the  Loch  Linnhe  trials.  This  suggests  that  the  presence  of 
the  mixed  layer  probably  renders  the  wake  undetectable  by  SAR  using  current 
technology. 

If  the  mixed  layer  thickness  is  increased  by  one  additional  meter,  the  result  is  shown  in 
Figure  31.  The  surface  velocity  amplitudes  are  reduced  by  a  factor  of  approximately  two. 
This  demonstrates  that  mixed  layers  above  the  main  internal  layer  have  a  very  negative 
effect  on  wake  detectability  at  the  surface.  The  wake  angle  is  increased  to  about  28°. 
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Figure  30.  The  velocity  profile  showing  the  effect  of  a  4  m  thick  mixed  layer. 
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Figure  31.  The  velocity  profile  showing  the  effect  of  a  5  m  thick  mixed  layer. 


An  example  of  a  simulation  main  window  with  the  Gibraltar  layer  profile  and  a  Wigley 
hull  representing  the  Sir  Tristram  moving  at  5  m/s  is  shown  in  Figure  32.  The  wake 
comprises  a  ‘V’  with  virtually  no  visible  internal  structure.  The  maximum  surface 
velocity  amplitude  along  the  parallel  track  is  a  fraction  of  a  millimeter  per  second.  As  the 
speed  of  the  ship  increases,  the  amplitude  increases  but  quite  slowly,  so  that  the  wake 
might  not  be  detectable  even  at  much  higher  speeds.  Details  of  the  velocity  along  the 
parallel  track  are  shown  in  Figure  33. 
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b?  Diffuse  Internal  Layer  Simulation 
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Figure  32.  Simulation  main  window  with  Gibraltar  profile  and  standard  parameters. 
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Figure  33.  The  velocity  profile  using  the  Gibraltar  layer  profile. 

Discussion 

Given  the  uncertainties  in  the  flow  velocities  derived  from  the  observations  during  the 
Loch  Linnhe  trials,  mutual  comparisons  between  observations  from  the  trials  and  the 
simulations  suggest  that  the  models  provide  reasonable  predictions.  The  diffuse  model  is 
clearly  superior  to  the  discrete  model  from  both  a  theoretical  point  of  view,  since  it  more 
closely  represents  real  layer  properties,  and  from  a  practical  point  of  view  in  that  it 
provides  a  better  match  to  the  observations.  It  is  not  restricted  to  narrow  layers  of  the 
Loch  Linnhe  type. 
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Matching  of  sets  of  wakes  includes  their  crest  patterns  and  their  amplitudes.  The  crest 
patterns  are  an  excellent  match  as  has  been  demonstrated  before  [6]  as  well  as  in  the 
present  study.  The  matching  of  SAR  wake  images  and  velocity  patterns  is  more  difficult 
because  it  requires  a  connection  between  radar  scattering  coefficients  of  the  ocean  surface 
and  the  surface  flows.  This  has  been  bypassed  in  part  by  using  the  derived  surface  flows 
from  the  trials.  If  we  assume  that  there  is  an  approximately  linear  relation  between  the 
spatial  derivatives  of  the  surface  flow  velocity  and  the  scattering  coefficient,  as  in  [9], 
then  there  is  also  a  good  match  between  the  general  appearances  of  the  observed  radar 
wake  and  the  simulated  wake. 

The  wakes  from  several  hull  forms  have  been  compared.  It  has  been  found  that  the  shape 
of  the  hull  is  not  particularly  important  for  internal  wave  wakes.  On  theoretical  grounds 
the  hull  shape  should  not  affect  the  wake  crest  pattern  but  only  its  amplitudes.  However, 
the  effect  on  the  amplitudes  is  quite  minor;  this  is  consistent  with  findings  for  the  discrete 
layer. 

Changing  the  ship  dimensions  affects  the  amplitudes  within  the  wake;  these  are 
approximately  proportional  to  the  maximum  transverse  area  of  the  ship  (beam  times  draft 
times  a  mid-ships  coefficient).  The  length  of  the  ship  affects  how  the  waves  generated  by 
sources  and  sinks  interfere.  The  wake  amplitudes  are  not  a  strong  function  of  ship  speed 
As  the  speed  increases  the  wake  narrows  and  amplitudes  tend  to  grow  but  quite  slowly. 

It  is  important  to  note  that  the  simulation  estimates  of  the  amplitudes  very  close  to  the 
wake  edge  are  sometimes  not  reliable  but  this  does  not  affect  most  of  the  results  here 
because  the  amplitudes  inside  the  wake  are  of  the  same  order.  The  problem  is  related  to 
the  neglect  of  internal  shock  waves  and  possibly  other  non-linear  phenomena.  When  the 
Wigley  and  other  models  are  used,  the  source  distribution  is  spread  out  along  the  ship’s 
length  and  this  tends  to  broaden  any  shock  wave  at  the  wake  edge  and  avoids  large 
amplitudes  analogous  to  the  sonic  boom  of  a  super-sonic  aircraft.  However,  problems  can 
arise  when  the  single  source  or  doublet  are  employed. 

Wake  Detectability 

Observations  from  the  Loch  Linnhe  trials  suggest  that  internal  wave  wakes  can  be 
observed  in  SAR  imagery  when  the  surface  flow  velocities  are  several  centimeters  per 
second  and  the  wavelength  of  these  changes  is  of  the  order  of  10  m.  Using  current 
technology,  such  as  RADARSAT-2  and  Sentinel- 1,  the  lower  limit  for  wake  detectability 
is  probably  of  the  order  of  one  centimeter  per  second. 

If  the  wavelength  of  waves  is  of  the  order  of  100  m  or  more,  as  for  the  interior  of  the 
Gibraltar  wake,  the  internal  wave  wake  generated  on  a  this  type  of  layer  is  probably 
difficult  to  detect  by  SAR  with  current  technology.  This  is  even  in  the  absence  of  a  mixed 
layer.  The  Gibraltar  wake  mainly  comprises  a  ‘V’  shape  and  the  width  of  each  wake  arm 
could  be  small  if  a  ship  were  moving  at  high  speed.  Though  this  would  not  necessarily 
provide  a  large  increase  in  the  wake  surface  flow  velocities,  it  would  increase  their 
derivatives  and  this  could  make  the  wake  arms  detectable. 
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The  existence  of  a  mixed  layer  just  a  few  meters  deep  is  sufficient  to  reduce  the  surface 
flows  that  might  be  detectable  to  millimeters  per  second.  The  presence  of  a  mixed  layer  is 
very  common  because  winds  of  a  few  meters  per  second  will  create  non-linear  effects, 
such  as  wave  breaking  that  eventually  mix  the  water  in  the  upper  ocean  to  a  depth  of  the 
order  of  the  wavelength  of  the  surface  gravity  waves.  Therefore  it  is  readily  understood 
why  internal  wave  wakes  from  surface  ships  are  rare;  they  will  only  be  detectable  when 
conditions  are  conducive  to  the  creation  of  internal  layers  and  the  water  surface  has  been 
calm  for  a  considerable  time. 

Conclusions 

The  diffuse  model  of  an  internal  layer,  the  wake  simulation  theory  and  its  implementation 
in  computer  code  provide  results  that  are  consistent  with  Loch  Linnhe  trial  data  and 
previous  simulations  using  a  discrete  layer  model. 

It  has  been  demonstrated  that  the  presence  of  a  surface  mixed  layer  dramatically  reduces 
the  wake  surface  flows  when  the  depth  of  the  layer  is  greater  than  the  draft  of  the  ship  so 
that  the  ship  does  not  penetrate  the  internal  layer  itself. 

The  internal  wave  wake  is  not  sensitive  to  the  hull  shape  and  depends  primarily  on  the 
principal  dimensions  of  the  ship  and  its  velocity. 

An  internal  wave  wake  that  can  be  identified  as  such  by  its  internal  structure  can  only  be 
observed  when  special  conditions  are  met.  Thus,  a  strong  layer  must  exist,  it  must  be 
close  to  the  surface  and  the  surface  mixed  layer  must  be  thin.  Wakes  on  weak,  deep 
layers  could  be  observable  as  a  ‘V’  with  no  internal  structure.  Again  there  must  not  be  a 
significant  mixed  layer  at  the  surface.  The  likelihood  of  observing  a  wake  on  the 
Gibraltar  profile,  for  example,  is  enhanced  if  the  ship  is  traveling  at  high  speed. 

The  simulations  have  been  based  on  the  simplest  of  practical  models.  There  are  other 
effects  that  can  enhance  wake  detectability  by  increasing  the  surface  flows.  One  example 
is  the  reflection  of  natural  waves  from  a  ship  hull.  A  wave  carries  momentum  and  during 
the  reflection  process  momentum  can  be  transferred  from  the  ship  to  the  wave  and  greatly 
increase  the  wave  amplitude.  A  second  example  involves  a  current  shear,  which  can 
amplify  waves,  or  at  least  reduce  their  attenuation.  The  low  ship  velocity  regime  (below 
the  super-critical  limit)  should  be  examined. 

A  significant  omission  in  the  simulation  is  the  likelihood  of  an  internal  shock  wave  at  the 
wake  edges.  This  could  be  important  for  ships  generating  wakes  on  deep  layers,  such  as 
the  Gibraltar  layer  and  should  be  corrected. 

No  evidence  has  been  found  for  soliton  wakes  but,  if  they  occur,  their  appearance  would 
resemble  the  Gibraltar  profile  wake  with  two  arms  and  negligible  internal  structure.  There 
is  a  possibility  that  the  way  in  which  solitons  form  by  non-linear  coalescence  of  a  wave 
train  could  also  increase  the  wave  amplitudes.  This  is  an  interesting  theoretical  problem. 
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Appendix  1:  Sturm-Liouville  Theory 

In  the  present  context,  we  are  dealing  with  a  simplified  Taylor-Goldstein  equation,  which 
can  be  written  in  the  form: 

^%  +  Ap(z)Q  =  0.  (8) 

dz 

Here  the  parameter  X  is  the  square  of  the  horizontal  wave  number  and  the  function  p(z) 
involves  the  horizontal  wave  frequency  and  the  B-V  frequency.  In  the  standard  theory,  p 
is  always  positive  but  in  this  application  it  can  be  negative.  Therefore  it  is  useful  to 
determine  if  this  is  an  important  limitation.  The  problem  is  to  solve  this  equation  subject 
to  boundary  conditions  such  that  Q  is  zero  at  z  =  0  and  at  z  =  L.  This  leads  to  a  sequence 
of  eigenvalues,  Xn  and  eigenvectors  Q„- 


The  standard  approach  is  to  first  show  that  the  operator  d2ldz~  is  self-adjoint  or  Hermitian. 
This  is  easily  done  for  the  present  boundary  conditions.  In  general  we  need  to  show  that 
(here  the  solutions  are  real  so  that  the  complex  conjugates  are  actually  not  relevant): 


d2u 

dz2 


dz 


f  l 


Vo 


d2v 

dz2 


dz 


Integrating  the  left-hand  side  by  parts  we  have: 

L  r  ?  r  ~\L 


*  d  u  , 
v  — -dz 
dz2 


v 


du 

dz 


•  dv  *  du 
dz  dz 


dz . 


(9) 


(10) 


o  ^  ^  Jo  o 

The  first  tenn  on  the  right  is  zero  because  of  the  boundary  conditions.  Another  integration 
by  parts  shows  that  the  operator  is  indeed  self-adjoint. 


Secondly,  it  is  shown  that  the  eigenfunctions  are  orthogonal.  Consider  two 
eigenfunctions,  Qn  and  Qm.  Multiplying  the  Taylor-Goldstein  equation  by  Qm  and 
integrating  we  have: 

} Q.  =  -V  J ■  (ID 

o  az  0 

However,  reversing  the  roles  of  m  and  n,  we  also  have: 


~K\ P(z)QnQJz  ■ 

0 


(12) 


The  left-hand  sides  of  (11)  and  (12)  are  equal  by  virtue  of  the  self-adjoint  property  so 
that,  if  the  two  eigenvalues  are  different: 

L 

(K-K)SpQ„QJz  =  0.  (13) 

0 


Therefore  the  integral  must  be  zero  irrespective  of  the  sign  of  p.  This  implies  that  the 
eigenvectors  are  orthogonal,  if  orthogonality  is  defined  in  terms  of  the  integral,  which 
includes  the  weighting  function  p. 
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When  m  =  n,  (13)  is  satisfied  automatically.  Moreover,  except  under  very  special 
circumstances  (namely  when  m  =  n  and  the  integral  is  zero),  they  can  also  be  nonnalized 
and  rendered  orthononnal,  i.e. 

L 

\pQ„QA  =  S..-  (14) 

0 

In  practice  there  appear  to  be  no  problems  in  normalizing  the  eigenfunctions. 


34 


Appendix  2:  The  Crest  Pattern  and  Wake  Waves 

Figure  34  shows  the  geometry  of  wake  formation.  A  disturbance  at  O,  created  by  ship  S 
propagates  at  an  angle  6  to  a  position  in  the  wake  described  by  angle  [1. 


y/=0 

Figure  34.  Geometry  of  wake  formation  using  group  velocity  concept. 


In  the  far  wake  the  phase  is  given  by: 

y/  =  cor-ks  =  ( co-kcg)z 

=  k(c-cg)i 


(15) 


where  r  is  the  propagation  time.  Therefore,  if  the  crests  are  separated  by  2n  and  n  is  an 
integer  indicating  the  crest  number,  the  time  interval  is  given  by: 

2tui  . 

T  = -  .  (16) 

k(c-c  ) 


It  follows  that,  using  the  phase  condition  c  =  UcosO,  the  coordinates  in  the  wake  relative 
to  the  ship  (x,  y )  are  given  by: 

x  =  (U  —  c  cos  0)t 


2miU(\-ccg/U2) 
k(c  ~cg) 
v  =  c  rsin# 

^  S 

Imic  g^J\-  c1  /U2 
k(c-cg) 


(17) 
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The  phase  and  group  velocities  are  determined  from  the  dispersion  relation  and  each  of 
these  can  be  regarded  as  a  function  of  either  wave  frequency  or  wave  number.  Therefore 
the  positions  along  the  crests  can  be  plotted  out  directly  from  (17). 


Because  the  medium  is  dispersive,  the  wake  waves,  which  propagate  outwards  from  a 
disturbance  created  by  the  ship,  soon  become  close  to  sinusoidal  in  shape.  At  each  point 
within  the  wake  they  are  represented  by  their  amplitude  and  their  wave  vector.  The  wave 
vector  can  be  found  from  geometrical  considerations.  Therefore,  in  terms  of  the  angle  fi 
in  Figure  34,  we  have: 


tan  /?  =  — 

x 


1  ~ccg/U2  ‘ 


(18) 


Now  c  and  cg  are  functions  only  of  co  or  k,  so  that  this  can  be  inverted  numerically  to  give 
co  or  k  (which  are  related  by  the  phase  velocity).  The  direction  of  k  can  be  found  from  the 
Doppler  equation: 

c  =  U  cos  9  (19) 


Once  the  wave  vector  at  a  point  has  been  detennined,  the  phase,  if/,  can  be  found  from  the 
geometry: 

if/  =  -k.r  =  -kr  cos (j8  +  9) .  (20) 

However,  to  evaluate  the  double  integral  in  Appendix  3,  we  need  to  know  how  the  phase 
varies  with  9.  We  assume  that  the  distance  r  is  large  so  as  to  be  able  to  use  the  stationary 
phase  approach,  which  is  consistent  with  the  concept  of  group  velocity. 
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Appendix  3:  Integral  Evaluation 

The  vertical  component  of  the  velocity  is  given  by: 

q z(r')  = 


1  1 

r  Qo(z)  c 

fdF0^ 

(2^)3  J 

'k2-k2  0 

V  dz  J 

t/k . 


The  horizontal  component,  qr,  is  found  from  the  divergence,  which  is  zero: 

dqz  dqr 


■  +  ■ 


=  0. 


dz  dr' 

At  the  surface  qz  =  0  but  the  derivative  is  finite.  From  (21),  we  have: 

i  dq: 


q,’  = 


k  dz 


The  direction  of  the  horizontal  flow  is  given  by  the  angle  6. 
We  can  evaluate  qr : 
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k0.r-w/4) 


(21) 

(22) 

(23) 


(24) 


To  find  the  second  derivative  of  the  phase  we  differentiate  (a  subscript  indicates  a 
derivative): 

y/g  =  ker  cos(/?  +  0)-  kr  sin(/?  +  0) 
y/ge  =  kder  cos(/?  +  6))  -  2 ker  sin(/?  +  0)  -  kr  cos(/?  +  0) 

Noting  that  y/g  =  0  at  a  saddle  point,  the  sine  term  can  be  eliminated  using  the  first  of 
these  equations  and  we  have: 


(25) 
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The  derivatives  of  k  can  be  found  from  the  equation  c=U  cos  9.  Then  we  have: 

—  ke  =  -U  sin  0 
dk 


d~c  ,  2 


+  —  k 
dk  dk 


-U  cos  0 


Therefore,  noting  that 
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dk 
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we  have: 


(26) 
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(28) 
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ke  = 


kU  sin# 

c~cg 

kU  cos#  k3U2  sin2  #  d2c 
■ - +  ■ 
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Combining  these  results  yields: 

Wee  = 


r U cos#  kU 2  sin2  #  dc 
- +  - 
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(c-c  )  dk 


W  ■ 


(30) 
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Appendix  4:  Matrix  Approach 

The  Runge-Kutta  and  similar  methods  of  integration  can  introduce  small  errors.  These 
can  be  very  important  when  finding  the  eigenvalues  and  eigenvectors.  The  problem 
occurs  when  the  frequency  approaches  the  maximum  B-V  frequency  when  the 
eigenvectors  fall  to  extremely  small  values  at  both  ends  of  their  range;  this  leads  to  rapid 
variations  that  imply  that  the  integration  step  size  should  be  reduced.  The  results  tend  to 
contain  errors  that  make  it  impossible  for  the  search  for  an  eigenvalue  and  particularly  the 
eigenvector  to  succeed.  The  eigenvalue  may  be  reasonable  but  the  apparent  eigenvector 
may  contain  a  large  positive  or  negative  spike  at  the  end  corresponding  to  the  ocean  floor. 


The  R-K  algorithm  is  not  wholly  responsible  for  this  problem  but  uncertainties  over  its 
accuracy  can  be  avoided  by  turning  to  a  matrix  approach  in  which  the  layer  is  exactly 
treated  in  a  step-wise  manner.  Within  each  step,  the  Taylor-Goldstein  equation  can  be 
solved  explicitly  and  the  boundary  conditions  matched.  The  solution  at  the  nth  step  takes 
one  of  the  following  forms,  depending  on  the  sign  of  the  final  term  in  the  Taylor- 
Goldstein  equation,  specifically  the  ratio  of  the  BV  frequency  at  a  given  depth  to  the 
frequency: 

\Aek-z  +Be~Kz,  k;  <  0 

y  =  \  n  .  (31) 

[Acos{knz)  + Bsm{knz),  k2n  >  0 


Here  A  and  B  are  constants,  z  is  the  distance  from  the  start  of  a  step  and,  from  the  Taylor- 
Goldstein  equation,  kn  is  given  by: 


k:  =k- 


N: 


co 


•1 


(32) 


In  terms  of  the  amplitude,  y,  and  its  derivative  with  respect  to  depth,  y\  in  the  first  case 
we  have: 


yn  =a  +  b 

y'n  =  k„(A-  B) 


(33) 


Solving  we  have: 
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(10) 


Therefore,  if  the  distance  between  steps  is  d,  the  amplitude  and  its  derivative  at  step  n+ 1 
are: 
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These  equations  can  be  put  into  matrix  form: 
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(35) 


>„+13 
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V 

cosh(A:nd)  sinh (knd)/kn 
kn  sinh (knd)  cosh (knd) 


ry, ' 


UJ 


9 

The  result  for  the  oscillatory  solution  is  identical  when  we  take  kn  to  be  negative,  or  in 
other  words: 


' yn+x 1 

^cos (knd)  sin (knd)lkn  N 

( y  \ 

y  n 

K.y'n+1  J 

K-knsm{knd)  cos  (knd)^ 

Jn; 

(36) 


where  k„  is  now  real. 


To  determine  the  eigenvalues,  k,  and  the  corresponding  eigenvectors  a  search  can  be 
implemented  to  find  the  eigenvalue  that  matches  the  boundary  conditions  both  at  the 
surface  and  at  the  bottom.  The  eigenvector  is  then  the  sequence  of  amplitudes  at  each 
step. 
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Appendix  5:  Comparison  of  Discrete  and  Diffuse  Layer  Theories 

It  was  considered  to  be  useful  to  compare  the  theories  for  the  discrete  and  diffuse  layers 
to  ensure  that  they  were  consistent  with  each  other;  therefore  they  were  less  likely  to  be 
invalid.  The  discrete  layer  theory  was  based  on  the  velocity  potential,  which  provides  the 
surface  velocity  components  quite  directly,  whereas  the  diffuse  theory  relies  on  the  zero 
divergence  of  the  velocity  at  the  surface. 


In  the  case  of  a  discrete  layer  we  work  in  terms  of  the  velocity  potential.  Near  the  surface, 
for  both  the  discrete  and  diffuse  cases  we  expect  this  to  be  given  approximately  by: 

0  =  Fcosh(kz)e"ikr  ,  (37) 

where  F  is  a  function  of  the  other  parameters  including  the  dispersion  relation.  It  follows 
that  the  vertical  and  the  surface  velocities  are  given  by: 


v.  = - =  kF  sinh (kz)e  ,kr 

dz 

vr  =  <^~  =  ikF  cosh  (kz)e~,k'r 
dr 

It  follows  that  the  vertical  fluid  velocity  is  zero  at  the  surface  (z.  = 
components  are  finite. 


(38) 


0)  but  the  horizontal 


For  the  diffuse  layer,  the  solution  near  the  surface  is  proportional  to  dvz/dz.  For  either 
layer,  this  is  given  by: 

&  =  k2F cosh (kz)e~ik-r .  (39) 

dz 

This  is  finite  at  the  surface.  Now  comparing  the  calculation  with  that  using  the  divergence 
of  the  velocity  (as  for  the  diffuse  layer)  yields: 

vr  =  —  =  ikF  cosh (kz)e~,k'r .  (40) 

k  dz 

This  is  the  same  as  in  (38).  Therefore  the  two  approaches  are  consistent.  The  factor  of  k 
in  (40)  ensures  that  the  amplitude  of  the  horizontal  velocity  near  the  wake  edge  does  not 
go  to  infinity. 
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